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ABSTRACT 

In this paper we make further numerical experiments assessing an accuracy degener- 
acy phenomena reported by A. Rogerson and E. Meiburg [7]. We also propose a modified 
ENO scheme, which recovers the correct order of accuracy for all the test problems with 
smooth initial conditions and gives comparable results with the original ENO schemes for 
discontinuous problems. 
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and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. Research also supported 
by NSF Grant DMS88-10150 and AFOSR Grant 90-0093. 



1 Introduction 


ENO (essentially non-oscillatory) schemes were first developed by Harten and Osher [2], 
Harten, Engquist, Osher and Chakravarthy [3], [4], to solve a hyperbolic conservation law 

f U t + f{ u )x = 0 /II') 

( u(x, 0) = u°(x) ' ' 


whose solution may be discontinuous even if the initial condition u°(x) is smooth. The 
philosophy of ENO schemes is to use an adaptive stencil, based on local smoothness, to 
automatically avoid interpolations across discontinuities. As a result, a formally uniform 
high order of accuracy, measured by local truncation errors, and sharp, essentially non- 
oscillatory shock transitions can be obtained. Shu and Osher [8], [9] later proposed efficient 
ways of implementing ENO schemes based on numerical fluxes and a special class of TVD 
(total- variation-diminishing) Runge-Kutta type high order time discretizations, which uses 
a conservative flux difference 


1 

Ax 



( 1 . 2 ) 


to approximate f(u) x in (1-1) to high order, where the numerical fluxes f are evaluated using 
high order interpolation polynomials constructed from adaptive stencils. 

Even if the numerical examples in [3], [8], [9] are very impressive, a convergence theory for 
ENO schemes is still not available, due to the nonlinearity of the scheme. A more difficult 
problem is to prove any convergence rate in smooth regions, since total variation stability, 
a tool used very often for proving convergence of nonlinear schemes in one space dimension, 
does not provide such a rate of convergence. 

Recently, A. Rogerson and E. Meiburg reported in [7] that for a linear, constant- coefficient 


version of (1.1): 


u t + u x = 0 
u(x, 0) = u°(x) 


(1.3) 


the convergence rate is dependent upon the initial condition u°(x) (assumed smooth) and 
may be inferior to what is predicted by the local truncation error analysis. Furthermore the 
error reduction during mesh refinements is not uniform and there are cases where a refined 
mesh gives a larger error. The source of this phenomena is traced to the fact that a linearly 
unstable stencil is initially chosen by ENO in a large portion of the domain. It was known 
to Harten et al [3] that the use of linearly unstable stencils will lead to frequent stencil 
switchings and loss of one order of accuracy due to the failure of error cancellation in the 
conservative form (1.2) (they verified this by using the initial condition u°(x) — e~ x in (1.3) 
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so that the initial stencil chosen by ENO is always the linearly unstable downwind stencil). 
The results in [7] indicates that the accuracy degeneracy phenomena can be more serious. 

In this paper, we perform further numerical experiments, including the use of different 
ENO schemes (cell- averaged version [3], [4] and the flux-pointvalue verion [8], [9]) and different 
time discretizations (TVD Runge-Kutta [8], [9] and exact time evolution [4]), to assess the 
accuracy degeneracy phenomena more thoroughly. We also propose a modified ENO scheme, 
which (1) introduces no additional computational cost over the original ENO scheme; (2) 
recovers the correct order of accuracy for (1.3) with different smooth initial conditions; and 
(3) gives comparable results with the original ENO schemes for discontinuous problems. 

All the computations in this paper are performed on a Cray-2 supercomputer. The largest 
number of grid points used is 5120. 


2 Numerical Experiments with ENO Schemes 

We use the third order ENO scheme based on numerical fluxes, and the third order TVD 
Runge-Kutta time discretization, described in detail in [8], [9]. CFL number is chosen as 
^ = 0.6 and the terminal time in (1.3) is chosen as t — 4 (two time periods). As in [7], We 
choose the initial conditions u°(x) in (1.3) as periodic functions with period 2, in order to 
exclude possible degeneracy of accuracy due to boundary conditions. 

First we briefly revisit the two examples used in [7], i.e. u°(x) = sin(nx) and u°(x) = 
sin 4 ( 7 rx), using the more commonly used third order in space and time ENO scheme with 
CFL = 0.6. As in [7], we can see that ENO is as accurate as the linearly stable centered 
scheme for u°(x) = sin( irx) (Figure la), but for u°(x) = sin 4 (wx) the Li error of the ENO 
scheme no longer monotonically decreases with the mesh size, and the final asymptotic order 
of accuracy is around 1.97, which is less than 3, predicted by local truncation error analysis 
(Figure 2a). By monitoring the stencil as in [7], we can see that, for tt°(x) = sin(-Kx), the 
initial stencil chosen by ENO falls completely inside the linear stability range 0 and —1 
(Figure lb), and it changes very little after two time periods (Figure lc). On the other 
hand, for u°(x) = sin 4 ( 7rx), ENO initially picks the linearly unstable stencils 1 and —2 in 
a large portion of the domain (Figure 2b) and then switches its stencil violently in order to 
balance the linear instability effects (Figures 2c and 2d). To see the different behavior of 
ENO scheme and the stable and unstable linear schemes, in Figure 2e we plot the logorithm 
of the Li error versus time for 160 grid points (which is the worst case for ENO from Figure 
2a). We can see that for a short time the L x error of ENO grows exponentially, similar 
to that of the linearly unstable scheme, but later ENO quickly stablizes itself, apparently 
through the nonlinear effect of stencil switching, while the linearly unstable scheme simply 
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blows up. However, what is lost to the accuracy during the brief exponential growth is not 
recovered, resulting in a L\ error for ENO which for this case is almost two magnitudes 
larger than that of the linearly stable centered scheme. 

Our next example is u°(x) = e” x . This is a discontinuous initial condition since we 
enforce periodicity. We compute the L\ error in the smooth region —0.5 < x < 0.5 at 
t = 4. From Figures 3a and 3e we can see that the linearly stable centered scheme now 
performs very poorly, with large Gibbs oscillations near the discontinuity and only first 
order accuracy even in the smooth region, due to the pollution of Gibbs oscillations. For 
this particular example the pollution can be avoided by using an exact time evolution, but for 
general nonlinear problems this pollution is very typical to high order linearly stable schemes 
applied to discontinuous problems. ENO scheme behaves much better than the linear scheme, 
cf Figures 3a and 3f, however we still observe the non-uniform decaying of errors with the 
mesh refinements, similar to the previous case u°{x) — sin 4 (7rx). This example was used 
in [3], without periodicity assumption, to illustrate the self- correcting mechanism of ENO 
schemes. Exact time evolution was used in [3], which explains the difference there and in 
Figure 3a. ^From Figure 3b we can see, as expected, that the initial stencil chosen by ENO 
is almost always the linearly unstable 1, and from Figures 3c and 3d we again observe the 
violent stencil switching at later times. 

Other initial conditions, such as u°(x) = sm 10 (7rx), u°(x) = (x 2 — l) 4 enforcing period- 
icity, ii°(x) — e $in ^ x \ u°(x ) = sm(cos(7rx)), etc., have been tested as well. We observe this 
degeneracy phenomena, to different extent, for all those cases. 

One naturally suspects whether this phenomena is related to the Runge-Kutta time 
stepping. As a matter of fact the result for ENO does become better with an exact time 
evolution [3], especially for u°(x) = e~ x . However, for most of the cases we have tested, 
this difference is only quantitative: the accuracy degeneracy phemonena is still present. 
Figure 4a shows the result with exact time evolution for u°(x) = sin 4 ( 7rx). Comparing 
with Figure 2a, we can see that the error of ENO now decays monotonically with mesh 
refinements, and the final asymptotic order of convergence is around 2.5, which is better than 
the Runge-Kutta case but still shows the degeneracy of accuracy. This suggests that the 
root of this degeneracy is the ENO spactial operator, although Runge-Kutta time stepping 
may compound the problem. 

We have also tested the original cell-averaged ENO schemes in [3], with the same initial 
conditions. We have not observed essential differences. Compare, e.g., Figure 4b with Figure 
4a. For the equation (1.3), the only difference between the two versions of ENO schemes is 
the initial condition (cell- averages versus point values). 

Different spatial orders for ENO schemes are also tested: up to seventh order with exact 
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time evolution. The accuracy degeneracy phenomena is present for all orders, but is getting 
less serious when the order increases. Compare, e.g., Figure 4c with Figure 4a. At seventh 
order (Figure 4c) the ENO result is quite satisfactory in the practical range of the number 
of grid points. This suggests that for a fixed cost, a higher order ENO with fewer grid points 
may be preferred to a lower order one with more points. 

Finally, the effect of this phenomena with different terminal time is tested. Figure 4d 
(t = 8) should be compared with Figure 2a (t = 4). We clearly see qualitatively similar 
results. 


3 A Modification of ENO Schemes 


The philosophy of ENO schemes is to adaptively choose stencils so that interpolations across 
discontinuities can be avoided. In smooth regions there seems to be conceptionally no need to 
choose linearly unstable stencils in the stencil-choosing procedure. The numerical evidence in 
[7] and in the previous section suggests that it is even harmful to the accuracy. Consequently 
some strategy to use a simpler fixed stencil linearly stable scheme in the smooth regions and 
to use ENO near discontinuities seems to be appealing. Earlier work for such strategies has 
already been documented in e.g., [5] and [10], with a motivation to save computational costs. 
These approaches will also solve the accuracy degeneracy problems for the examples with 
smooth initial conditions in [7] and in the previous section, since a linearly stable centered 
scheme is used throughout the region except for possibly a few isolated cells near critical 
points, which seems not to affect the accuracy in the numerical tests (S. Osher, private 
communication). 

In this section we discuss a modification of ENO schemes which involves only a slight 
change in the coding, without increasing the computational costs. A similar modification 
has been used in [1] for a different purpose (getting smoother steady state solutions). Instead 
of using 


i(j) = j ( upwinding ) 

for k = 2, ..., r + 1 : 

if ( abs(H[i(j),k]).gt.aba(H[i(j ) - l,fc])) i(j) = i(j ) - 1 
end for 


(3.1) 


to determine the left-most point i(j) in the stencil for computing the numerical flux /•, 1 
in (1.2) (H[i, fc] is the k-th divided (or undivided during the actual coding) difference of the 
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function H(x), whose definition and relations to /(u(x)) can be found in [9]), we use instead 

i(j) = j (unwinding) 
for k = 2, ..., r + 1 : 

if (i(j).gt.ic(j)) then 

if (2 * abs(H[i(j),k]).gt.abs(H[i(j) - l,As])) *(j) = i(j) ~ 1 ^3 2 ) 

else 

if (abs(H[i(j),k]).gt. 2 * abs(H[i(j) - 1, A:])) i(j) = i(j ) - 1 
end if 
k end for 

where i c (j ) is the left-most point in the linearly stable centered stencil. By using this 
modification, we are trying to bias towards the linearly stable centered stencil as much as 
possible in the smooth regions. Since the ratio of the two differences being compared in (3.1) 
is 1 + 0( Ax) in smooth regions, it is easy to prove that in smooth regions where all the 
derivatives of H(x) are non-zeroes, (3.2) will give the centered stencil when Ax is sufficiently 
small. Near shocks the modifications should not affect the non-oscillatory philosophy, since 
the differences being compared are then of different magnitudes. 

We first test the performance of the modified ENO scheme (3.2) on the examples in [7] 
and in the previous section. As one can see from Figures la and 2a, the modified ENO is 
fully third order accurate and is in fact indistinguishable from the linearly stable centered 
scheme for smooth initial conditions. From Figures 3a and 3g we can see that the modified 
ENO performs much better than the linearly stable centered scheme and better than the 
ENO scheme, by producing fully third order accuracy in smooth regions and a non-oscillatory 
shock transition, when the initial condition is discontinuous. From Figures 2f, 2g and 3h, 3i 
we can see that, in the smooth regions, the stencils of the modified ENO scheme are almost 
always within the linearly stable choices 0 and —1. This is probably the main reason for its 
excellent behavior. We have also tested the modified ENO scheme in many other cases, such 
as all those mentioned in the previous section. Full high order of accuracy predicted by local 
truncation error analysis is always observed. 

We remark here that (3.2) with r = 1 is a MUSCL type second order (in the L\ sense) 
TVD scheme discussed in [6] (with the minmod function replaced by a minimum- in- absolute- 
value function). This is a major reason for our choice of the factor 2. Any factor bigger than 
1 should have the same effect of biasing towards a centered stencil asymptotically. We do 
not recommend a factor bigger than 2 since it will not produce a TVD scheme for r = 1. 

A more important test for the performance of the modified ENO scheme is to apply 
it to problems with both discontinuities and detailed smooth structures, we have already 
seen the application to (1.3) with u 0 (x) = e~ x in Figures 3a and 3g. We now apply it 
to the one dimensional Euler equations of gas dynamics, which is a nonlinear hyperbolic 
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system. We choose two shock tube problems and a problem of shock interaction with density 
perturbations, see [9] for details of equations and parameters. Comparisons between Figures 
3f and 3g, 5a and 5b, and 6a and 6b, reveal almost identical performance of ENO and 
modified ENO schemes. For the more interesting test case of shock interaction with density 
perturbations, Figures 7a and 7b, we observe a slightly better resolution of the modified 
ENO versus ENO with 160 grid points, which is margined for ENO to resolve this structure. 
However, the modified ENO has a slight over- compression effect (step phenomena) in the 
smooth regions, which is more visible in Figure 7d with 400 grid points. This is anticipated 
since (3.2), in many cases (e.g. for r = 1), is similar to the artificial compression used in [10] 
and [9]. 


4 Concluding Remarks 

ENO schemes may lose the full high order accuracy predicted by local truncation error 
analysis, for scalar linear conservation laws with smooth initial conditions. A modified 
ENO scheme, which does not increase the computational cost, can overcome this accuracy 
degeneracy problem for the test problems. For systems of nonlinear conservation laws with 
solutions containing both shocks and detailed smooth structures, ENO and modified ENO 
perform comparably, with the modified ENO having a slightly better resolution but showing 
some over- compression effects. 

Acknowledgment: We thank Stanley Osher for many valuable suggestions. We also thank 
David Gottlieb, Eckart Meiburg and Andrey Rogerson for helpful discussions. 
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NUMBER OF GRIO POINTS 


Figure la: L x error as a function of the number of grid points. u°(x) — sinfax), t = 4, 
££ = 0.6. Solid line: linear centered; Dashed line: ENO; Dot-Dashed line: modified ENO. 
Number of grid points: 10,20,40,80,160,320,640,1280,2560,5120. Log-log scale: slope 
represents order of accuracy 
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Figure lb: The stencil choice of ENO as a function of the grid index. u°(x) = sin(rrx), 
t = 0. The possible stencils for computing the flux / J+ i are:— 2 = 

— 1 = (xj_ 2 , Xj+i); 0 = ( Xj-i,Xj,Xj + i,Xj+ 2 ); 1 = ( Xj,Xj+i,Xj + 2 ,Xj +3 ); 0 (centered) 
and -1 (upwinding) are hnearly stable stencils, 1 and -2 are linearly unstable stencils 
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Figure lc: same as lb except that t = 4 
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Figure 2e: L x error, in logrithm scale, as a function of time t, for u°(x) = sm 4 ( 7 rx), with 
160 grid points. Solid line: linear centered; Dashed line: ENO; Dot-Dashed line: modified 
ENO; Dotted line: linearly unstable scheme 
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Figure 2f: same as 2b except that the modified ENO is used 
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Figure 3d: same as 3b except that t = 4 



Figure 3e: t — 4, u°(x) = e -x . Solid line: exact solution; Circles: solution of the linear 
centered scheme 
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Figure 4c: same as 4a except that seventh order in space is used 
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NUMBER OF GRID POINTS 


Figure 4d: same as 2a except that t — 8 





Figure 5a: shock tube problem with (pl,ul, Pl)=(1j 0, l)i [PRt u R> Pr)=(0.125, 0,0.1). Den- 
sity is plotted. 100 grid points. Solid line: exact solution; Circles: ENO 



Figure 5b: same as 5a except that Circles: modified ENO 
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Figure 6a: same as 5a except that ( Pl> u l , Pl)=(0.445, 0.698, 3.528), ( Pr,ur , Pr)— (0.5, 0, 0.571) 



Figure 6b: same as 6a except that Circles: modified ENO 
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Figure 7 a: shock interaction with density perturbations. p L )=(3.857143, 2.629369, 10.333333), 

{pRy^Ry + 0.2sin(5x), 0, 1). Density is plotted. Solid line: a converged solution com- 

puted by ENO with 1200 grid points; Circles: ENO with 160 grid points 



Figure 7b: same as 7a except that Circles: modified ENO with 160 grid points 
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